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GENERAL  NON-ITERATIVE  SOLUTION  OF  THE 


INVERSE  AND  DIRECT  GEODETIC  PROBLEMS 


I.  INTRODUCTORY  BACKGROUND 

Earlier,  at  the  Array  Map  Service,  the  writer  published  a  com¬ 
prehensive  study  [1]  for  the  rigorous  non-iterative  inverse  solution 
of  long  geodesics,  following  its  presentation  at  the  XI  General 
Assembly  of  the  International  Union  of  Geodesy  and  Geophysics.  (At 
present,  a  translation  by  the  German  Geodetic  Commission  and  an 
abstract  by  the  U.  S.  S.  R.  Academy  of  Sciences,  complete  with  final 
formulas,  are  also  available.)  The  method  received  favorable  commen¬ 
taries  from  several  authoritative  writers  on  the  subject.  A  com¬ 
parative  evaluation  study  published  in  [2]  indicated  that  the  method 
was  as  accurate  as  contemporary  solutions,  yet  simpler  and  shorter 
to  compute.  The  procedure  does  not  require  any  special  geodetic 
tables  and  calls  for  relatively  few  trigonometric  interpolations. 

Of  greatest  interest  was  the  fact  that  it  was  the  first  rigorous 
non-iterative  inverse  solution  to  go  not  only  beyond  the  first  power 
of  spheroidal  flattening  but  through  all  the  cubic  terms  of  flatten¬ 
ing.  Yet,  it  was  practical. 

II.  THEORETICAL  REFINEMENT  OF  NON-ITERATIVE  INVERSE 

Since  [1]  represented  only  the  formative  stage  of  development, 
the  next  phase  consisted  of  making  a  thorough  mathematical  analysis 
of  the  formulas  in  order  to  determine  whether  there  were  any  con¬ 
cealed  intrinsic  properties  or  relationships  which  could  be  used  to 
obtain  an  optimum  solution.  This  analysis  resulted  in  the  uncover¬ 
ing  of  three  basic  quantities,  a,  m,  and  0,  the  exclusive  substi¬ 
tution  of  which  promised  a  concise,  orderly  pattern  for  the  terns 
of  the  two  main  spheroidal  power  series,  x  and  S,  that  were  origi¬ 
nally  given  in  [1],  at  the  top  of  page  18  and  the  bottom  of  page 
19,  respectively. 

III.  DEVELOPMENT  OF  CORRESPONDING  DIRECT 


Since  there  were  indications  that  a,  m,  and  0  were  a  rather 
unique  set  of  quantities  (as  will  be  shorn  later),  an  attempt  was 
also  made  to  introduce  their  equivalent  into  the  formulation  of  a 
corresponding  Direct  solution,  which  so  far  was  lacking.  As  if  by 
design,  two  spheroidal  power  series  resulted,  again  with  a  simple, 
orderly  pattern  of  terms.  Moreover,  the  format  was  identical  to 
that  of  the  Inverse.  The  three  corresponding  basic  quantities  were 
denoted  as  a1  ,  ,  and  '0  .  This  Direct  solution  was  derived  from 

the  formulas  on  pages  14^and  15  of  [1],  Essentially,  the  power 
series  of  the  quantity  S  (now  appearing  also  in  Appendix  E)  was  used 


Co  solve  for  0(J ;  Chen  Che  expression  ol  0Q  was  in  Curn  placed  inCo 
Che  \  power  series.  A  critical  factor  in  Che  mathemacical  determi¬ 
nation  of  a,  and  it^  was  the  proper  choice  of  their  common  smaller 
order  term  .  5e'3  sin2  ^  ,  which  is  probably  one  in  a  series  of  others 
required  for  the  orderly  theoretical  extension  of  the  solution  to 
higher  degree. 


TV.  NOTES  ON  COMPUTATION  FORMS 

The  above  efforts  led  to  the  Inverse  and  Direct  computation  forms 
now  shown  in  Appendices  A  and  B,  respectively.  The  main  Inverse 
spheroidal  expressions  are  indicated  by  (S  -f  bQ)  and  (\  -  L)  -f  c, 
while  the  Direct  by  0e  and  (L  -  \)  -f  cos  g0 .  It  is  to  be  noted  that 
their  outer  coefficients  are  simple  product  combinations  of  a  and  m 
or  &x  and  ,  while  their  bracketed  expressions  are  functions  of  only 
the  variable  0  or  0^  and  powers  of  the  spheroidal  parameter  f  or  e'2. 
By  tabulation  of  the  functions  of  0  and  0.g  so  that  they  may  be  rapid¬ 
ly  obtained  by  interpolation,  it  would  then  be  very  simple  to  multi¬ 
ply  them  by  the  easily  determined  outer  coefficients.  Since  con¬ 
ventional  trigonometric  tables  can  be  used  to  obtain  the  first-order 
term  0,  and  no  tables  are  required  to  obtain  the  first-order  term 
0  ,  only  short  tables  for  the  second-  and  third-order  terms  need  be 
d^awn  up.  The  third-order  terms,  which  are  small,  vary  sufficiently 
slowly  for  visual  interpolation. 

In  addition  to  their  two  pairs  of  main  power  series  and  the 
basic  quantities  necessary  to  numerically  evaluate  them,  these  In¬ 
verse  and  Direct  computation  forms  generally  contain  one  simple 
closed  trigonometric  expression  for  each  required  final  quantity  ex¬ 
cept  when,  for  example,  the  trigonometric  cofunction  is  given  as  an 
alternative  for  occasions  when  a  weak  determination  or  unlinear  in¬ 
terpolation  may  otherwise  result.  The  forms  may  appear  cluttered 
with  rules  for  choosing  signs,  trigonometric  quadrants,  and  so  forth. 
Actually,  these  do  not  entail  any  added  calculations  but  simply  de¬ 
fine  the  problem  without  ambiguity.  In  addition,  the  choices  pro¬ 
vide  for  greater  generality  and  more  varied  applications.  For  ex¬ 
ample,  one  may  calculate  either  the  shorter  geodesic  between  two 
given  positions,  or  the  longer  geodesic  around  the  spheroid's  back 
side.  Also,  the  subscripts  1  or  2  may  be  assigned  to  either  position 
without  fear  of  ambiguity.  For  less  accuracy,  terms  in  f2  and  e'4 
may  be  omitted. 

V.  FORMULAS  FOR  VERY  SHORT  AND  LONG  GEODESICS 

From  its  inception,  the  non-iterative  derivation  was  developed 
primarily  for  very  long  geodesics;  therefore,  the  auxiliary  trigo¬ 
nometric  functions  which  were  correspondingly  designed  provided  the 
greatest  simplicity  at  the  expease  of  generality.  Recently,  how¬ 
ever,  this  Agency  as  well  as  the  Amy  Map  Service  placed  a  requirement 


for  a  single  set  of  formulas  applicable  to  very  short  as  well  as  very 
long  lines,  since  they  were  to  be  used  also  in  the  adjustment  of  tri- 
angulation  and  trilateration  ground  nets.  It  was  felt  that  the  basic 
long  line  formulas  given  in  Appendices  A  and  B  appeared  particularly 
convenient  for  the  electronic  computers  which  were  to  make  the  calcu¬ 
lations.  However,  in  order  to  obtain  the  same  or  even  greater  accu¬ 
racy  for  very  short  geodesics,  the  alternate  formulas  presented  in 
Appendix  C  were  provided.  The  reason  for  the  increased  accuracy  re¬ 
quirement  is  easily  understood  when  one  considers,  for  example,  that 
if  the  length  of  a  line  is  decreased  a  thousandfold,  the  positions 
of  the  new  endpoints  must  be  known  a  thousand  times  more  accurately 
to  maintain  a  constant  azimuth  accuracy.  This  means  that  the  lati¬ 
tude  and  longitude  will  have  a  greater  number  of  decimals  and,  there¬ 
fore,  additional  significant  digits.  In  order  to  avoid  the  carrying 
of  too  many  fixed  places,  which  are  more  apt  to  be  affected  by  round¬ 
ing  errors  if  there  is  no  spare  digit  capacity,  Appendix  C  provides 
formulas  whose  terns  are  generally  very  small  when  a  geodesic  is 
very  short,  so  the  computations  can  be  done  conveniently  by  floating 
point  for  greater  decimal  accuracy.  To  obtain  the  full  required 
accuracy,  no  additional  terms  need  be  added  to  the  power  series 
formulas,  because  (as  shown  in  the  second  paragraph  of  Appendix  F) 
they  converge  to  more  good  decimal  places  for  shorter  geodesics. 

Since  many  of  the  small  quantities  in  Appendix  C  consist  of  sines  of 
small  angles,  their  evaluation  is  especially  adaptable  to  electronic 
computers,  which  by  means  of  floating  point  can  readily  calculate 
trigonometric  series  that  inherently  converge  to  additional  decimals 
for  such  small  angles.  Actually,  the  formulas  in  Appendix  C  are 
equally  applicable  to  short  and  long  lines,  so  only  one  set  of 
equations  need  be  programmed  into  an  electronic  computer.  A  floating 
point  formula  for  a  more  accurate  cosine  of  large  absolute  latitude 
is  also  included  in  Appendix  C. 

VI.  CONCLUDING  REMARKS 

Appendix  D  provides  the  complete  numerical  calculations  for  a 
very  short  and  a  very  long  geodesic--!  mile  and  6,000  miles,  re¬ 
spectively.  In  each  case,  the  Direct  solution  provides  a  check  on 
the  Inverse.  The  discrepancies  between  the  two  types  of  solutions 
are  given  at  the  end  of  Appendix  D.  The  better  positional  accuracy 
provided  for  the  short  line  by  the  formulas  of  Appendix  C  is  con¬ 
vincingly  shown. 

Appendix  E  provides  a  non-iterative  inverse  solution  of  higher 
order  accuracy  (that  is,  through  f3  and  e,s  terms)  for  use  as  a 
theoretical  check  on  Direct  or  other  Inverse  formulas.  In  Appendix 
F,  several  interesting  types  of  inter-relations  of  the  terms  of  the 
power  series  are  discussed  and  illustrated.  These  include  relation¬ 
ships  between  numerical  coefficients  as  well  as  algebraic  terms.  In 
Appendix  G,  meridional  arc  formulas  are  derived  as  special  cases  of 
the  Inverse  and  Direct. 
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In  conclusion,  it  should  be  noted  that  the  Inverse  case  of 
almost  antipodal  positions,  treated  on  pages  24  through  25  of  1], 
is  omitted  here  because  of  its  rare  practical  occurrence.  Also,  the 
elimination  of  3  by  substitution  in  terms  of  the  given  B  is  not 
undertaken  because  simple  closed  functions  herein  would  become  series 
expansions. 


M. » t  wf  M* vii  Uti  iO»  >|v^  IMIM  tnuMIM 


was*  JitobJ  tiaibi  -Sinrki 


BIBLIOGRAPHY 


1.  Sodano,  E.  M. :  "A  Rigorous  Non-Iterative  Procedure  for  Rapid 

Inverse  Solution  of  Very  Long  Geodesics," 
Bulletin  Geode's ique,  Issues  47/48,  1958, 
pp.  13  to  25. 


2.  Bodemiiller,  H. :  Beitrag  zur  Lflsung  der  zweiten  geodatischen 

Hauptaufgabe  fUr  lange  Linien  nach  den 
direkten  Verfahren  von  E.  Sodano  und  H. 

Moritz  sowie  nach  dem  sog.  Einschwenkverfahren 
Braunschweig,  im  Juni  1960. 


I 


5 


APPENDICES 


Appendix 

A 

B 

C 

D 


Item 

INVERSE  COMPUTATION  FORM 

DIRECT  COMPUTATION  FO «M 

ALTERNATE  INVERSE  AND  DIRECT  FORMULAS 

NUMERICAL  ILLUSTRATIONS  OF  INVERSE  AND 
DIRECT 


E  THEORETICAL  FORMULAS  FOR  HIGHER  ACCURACY 

F  INTER-  RELATIONS  OF  THE  TERMS  uF  Tilt 

POWER  SERIES 

G  MERIDIONAL  ARC  AS  SPECIAL  CASE  OF  NON¬ 

ITERATIVE  INVERSE  AND  DIRECT 


6 


APPENDIX  A 


INVERSE  COMPUTATION  FORM 


Given:  B1  ,  L:  =  Geodetic  latitude  and  longitude  of  any  point. 

Bg ,  Lg  =  Latitude  and  longitude  of  any  other  point. 

(South  latitudes  and  west  longitudes  considered  negative.) 
Required:  a,  S  =  Geodetic  azimuths  clockwise  from  north  and  distance. 


aQ ,  bQ  =  Semimajor  and  semiminor  axes  of  spheroid. 

bo 

f  =  Spheroidal  flattening  =  1  - 


L  =  (Lg  -  ^  )  or  (Lg  -  1^  )  +  [sign  opposite  of  (Lg  -  1^)]  (360°) 

Use  whichever  L  has  an  absolute  value  <  or  >  180°,  according 
to  whether  the  shorter  or  the  back-side's  longer  geodesic  is 
intended.  However,  for  meridional  arcs  ( I L '  =0°  or  180°  or 
360°),  use  either  L  but  consider  it  (+)  for  the  shorter  and 
(-)  for  the  longer. 

tan  3  -  (tan  B)  (1  -  f)  when  I  Bl  £  45° 
or  cot  g  -  (cot  B)  t  (1  -  f)  when  - Bl  >  45° 

I 

I 

a  =  sin  3X  sin  33  ;  b  =  cos  3^03  3S  ;  cos  0  =  a  +  b  cos  L. 

sin  0  =  -  y^sirTL  cos  3JP  +  (sin  3ffl  cos  31  -  sin  3,  cos  33  cos  Lp 

The  sin  0  is  (+)  for  the  shorter  arc  and  (-)  for  the  longer. 
Compute  the  radical  entirely  by  floating  decimals  to  prevent 
loss  of  digits,  especially  for  very  short  geodesics. 

0  =  Positive  radians  in  proper  quadrant,  reference  angle  being 
determined  from  sin  0  or  cos  0,  whichever  has  the  smaller 
absolute  value. 

c  -  (b  sin  L)  r  sin  0;  m  =  1  -  c3 . 
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=  [(1  +  f  +  fS)  0] 


-O 

+  a  [(f  +  f2)  sin  0  -  (~— )  03  esc  0] 

+  m  [-  (L±JL)  0  -  (1-i-L)  sin  0  cos  0  +  (|1)  02  cot  0] 

f2 

+  a3  [-  (L_)  sin  0  cos  0] 

+  m3  [(|^)  0  +  (~)  sin  0  cos  0  -  (~)  03  cot  0  -  (|_)  sin  0  cos3  0] 


-f  am  [(-— )  02  esc  0  +  (— )  sin  0  cos3  0] 

2  2 


^  -  [<£  +  £3>  0] 

f2 

+  a  [-  (— )  sin  0  -  (f3)  02  esc  0] 

5  f2  f2 

+  m  [-  (— — )  0  +  (-— )  sin  0  cos  0  +  (f2)  02  cot  0]  radians 

cot  ^  a  =  (sin  ps  cos  px  -  cos  \  sin  p:  cos  p3 )  t  sin  \  cos  p3 

cot  Og  x  =  (sin  Pg  cos  px  cos  \  -  sin  p1  cos  P3 )  t  sin  \  cos  g. 

For  meridional  arcs,  consider  a  as  having  0°  reference  angle,  and  ob¬ 

tain  only  the  signs  of  the  cotangents  by  disregarding  the  denomina¬ 
tors.  For  other  geodesics,  replace  cotangent  by  tangent  when 
| cot  aj  >1,  by  taking  the  reciprocal  of  the  quotient's  value. 


Quadrant  of  _3 

'Quadrant  of  a3  _x 

if  L  is  (+) 

...and  cot  (tan)  of  _3 

is  (+)  or  (-),  Q-^g  is  in 
quad  I  or  II,  respectively. 

...and  cot  (tan)  of 
is  (+)  or  (-) ,  o-g  _x  is  in 
quad  III  or  IV,  respec¬ 
tively. 

If  L  is  (-) 

...and  cot  (tan)  of  ax  _g 

is  (+)  or  (-),  q^.2  is  in 
quad  III  or  IV,  respec¬ 
tively. 

...and  cot  (tan)  of  a3 
is  (+)  or  (-),  oi2_x  is  in 
quad  I  or  II,  respec¬ 
tively. 
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APPENDIX  B 

DIRECT  COMPUTATION  FORM 

Given:  B1 ,  =  Geodetic  latitude  and  longitude  of  any  point  1. 

cvx  _ 2 ,  S  -  Azimuth  clockwise  from  north  and  distance  to 
any  point  2. 

Required:  Geodetic  p  Bj,  and  Lg . 

(South  latitudes  and  west  longitudes  considered  negative.) 


aQ ,  bc  =  Semimajor  and  semiminor  axes  of  spheroid. 

f  =  Spheroidal  flattening  =  1  -  — - 

aQ 

e'2  =  Second  eccentricity  squared  =  (a2  -  b")  -f  b3 

tan  0  =  (tan  B)  (1  -  f)  when  1 b!  £  45° 
or  cot  8  =  (cot  B)  t  (1  -  f)  when  I  b!  >  45° 
cos  80  =  cos  gx  sin  ax  _s  ;  g  =  cos  gx  cos  al_s', 

i%  =  (1  +  ~~  sin2  g1)  (1  -  cos2  80);  0g  =  (S  t  b0)  radians; 
,  2 

ax  =  (1  +  ---  sin3  px)  (sin2  0X  cos  0g  +  g  sin  px  sin  0g). 


0n  [0o] 


1 2 


+  ai  ["  ~  sin  0q] 


2  rS- 

la 

T  '  4 


+  "h  [-  V"  +  hr  sin  cos  0J 


?  r  5  6  *  ^ 

+  a3  [  4_~  sin  0  cos  0  ] 
0  b  b 


-  rile'4 

+  K  i-nr~  0 


64  64 


13e'4  e'4  a  5e 14  .  3  . 

^  -L - sin0  cos30  ] 


sin0gCos0g  -  ■  g  ’  0gCos  0g  -}~  ^ 


+  ax  mx  [ 


qe  14  e'4  5e 14  2  i  , 

_ —  sin  0  +  -  0  cos  0  - - sin  0  cos  0  radians 

8  S  4  S  S  8  S  SJ 
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c°t  ara _2  =  (g  cos  0O  -  sin  px  sin  0Q)  -:-  cos  pQ 

For  meridional  arcs,  consider  ag _x  as  having  0°  reference  angle, 
and  obtain  only  the  sign  of  the  cotangent  by  disregarding  the  de¬ 
nominator.  For  other  geodesics,  replace  cotangent  by  tangent 
when  I  cot  o-g  _x  I  >  1,  by  taking  the  reciprocal  of  the  quotient's 


value. 

Quadrant  of  g  1 

If  (0°£  C^.g  £  180°) 

....and  cot  (tan)  of  as _x  is  (+)  or  (-), 
as  x  is  in  quad  III  or  IV,  respectively. 

If  (180° <  _a  <  360°) 

_ and  cot  (tan)  of  a2 _x  is  (+)  or  (-), 

cts  is  in  quad  I  or  II,  respectively. 

cot  \  =  (cos  j3x  cos  0Q  -  sin  fa  sin  0O  cos  _2 )  r  sin  0Q  sin  ax  _3 

For  meridional  arcs,  consider  \  as  having  0°  reference  angle,  and 
obtain  only  the  sign  of  the  cotangent  by  disregarding  sin  g. 
For  other  geodesics,  replace  cotangent  by  tangent  when 
I  cot  \  I  >  1,  by  taking  the  reciprocal  of  the  quotient's  value. 


Quadrant  and  Sign  of  \ 

When  O°<0e^18O° 

(sin  0Q  considered 
positive) 

When  18O°<0o<;36O0 
(sin  0O  considered 
negative) 

and 

(o^.^iso0) 

...then  if  cot  (tan)  of  \ 
is  (+)  off  (-),  \  is  in 
quad  I  or  II,  respect¬ 
ively. 

...then  if  cot  (tan)  of  X 
is  (+)  or  (-),  X  is  in 
quad  III  or  IV,  respect¬ 
ively. 

and 

dsc^i^oeo0: 

...then  if  cot  (tan)  of  X 
is  (+)  or  (-),  the  associ¬ 
ated  angle  is  in  quad  III 
or  IV,  respectively,  and 

X  is  obtained  by  sub¬ 
tracting  360°. 

...then  if  cot  (tan)  of  X 
is  (-f)  or  (-),  the  associ 
ated  angle  is  in  quad  I 
or  II,  respectively,  and 

X  is  obtained  by  sub¬ 
tracting  360°. 

1 
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[-  f  0J 

b 

[— _  sin  0  ] 

L  2  WSJ 

f  — 0  -  3f2.  sin  0  cos  0„ ]  radians 

"  4  S  4  b  S 

Ifc  =  h  +  L 

[If  1 3ug  I  >180°,  modify  Lj  by  adding  _or  subtracting  360°,  according 
to  whether  it  is  initially  negative  or  positive.] 


sin  p2  =  sin  cos  0O  +  g  sin  0Q 


cos  p3  =  V(cos  30)2  +  (8  cos  0o  “  sin  Pi  sin  0o )2 

Compute  the  radical  entirely  by  floating  decimals  to  prevent 
loss  of  digits,  especially  for  large  absolute  latitudes. 

tan  ga  =  (sin  p2  4-  cos  pg ) 

Use  whichever  has  the 

or  cot  3a  =  (cos  ^2  4  sin  3a )  smaller  absolute  value. 

Obtain  tan  (or  cot)  of  b  from  earlier  defined  relation  of  B  to  3. 

2 

Determine  (-90°sil3  ^909 ) ,  applying  sign  of  its  tan  (or  cot). 


cos  p0 

+  a1 


APPENDIX  C 

ALTERNATE  INVERSE  AND  DIRECT  FORMULAS 
(For  very  short;  as  well  as  long  geodesics) 

The  following  alternate  formulas  for  corresponding  ones  in 
Appendices  A  and  B  are  designed  to  maintain  or  appropriately  in¬ 
crease  the  accuracy  of  various  elements  of  short  geodesics,  without 
decreasing  the  accuracy  of  long  geodesics.  The  formulas  specifically 
take  advantage  of  inherently  small  quantities  and  of  small  differ¬ 
ences  of  given  large  quantities,  so  as  to  provide--through  the 
application  of  floating  point  calculations--increased  decimal  place 
accuracy  without  requiring  additional  operational  digits.  The  small 
angles  involved  are  especially  adaptable  to  electronic  computers, 
which  by  means  of  floating  point  can  readily  obtain  greater  decimal 
accuracy  inherent  in  trigonometric  power  series  of  such  small  angles. 

1.  FOR  INVERSE  SOLUTION; 

sin  0  =  -  y(sin  L  cos  j3s)2  +  [sin(0a-  03 )  +  2  cos  fa  sin  $3  sin2  ^]2 

cot  a  -  [sin(3a-  )  +  2  cos  p2  sin  sin3  |]  ••  cos  02  sin  \ 

cot  asj.  =  [  s  in  ( 02  -  )  -  2  cos  0X  sin  fa  sin3  1]  *  cos  ^  sin  \ 

where 

(P8-Pl)  =  (Ba-B.)  +  2  [sin(Bs-B1)]  [(n  +  n3  +  n3 )  a  -  (n  -  n3  4-  n3)  b] 


2.  FOR  DIRECT  SOLUTION: 

b3  =  %  +  (0a -Pi)  +  2  [•3in(33-p1 )]  [(n  +  n3)  cos  (Bg+I^)  +  n3cos(B3  )  ] 


oIu  Vp3  ‘Pi  ) 
md  the  required 


—  sin  0£  cos  _2  -  2  siu3 
approximate  Ba  and  cos  ga 


2  siu  Pi  cos  *2 

are  obtained  in  Appendix  B. 


[*•. 
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3.  FOR  INVERSE  AND  DIRECT  AT  GIVEN  ABSOLUTE  LATITUDES  >  45° r 

cos  0  =  sin  ((90  +  B)  -  2  [sin  (90  +  B)]  (n  4-  n8  +  n3)  sin  3} 

the  upper  and  lower  signs  of  which  are  applied  for  the  northern  and 
southern  hemispheres,  respectively. 

In  the  preceding  three  sets  of  formulas,  n  =  (a(J  -  b0) 

(ao  +  t>0).  Some  smaller  coefficients  of  the  almost  negligible  n3 
have  been  removed  because  they  are  unsymmetric,  and  because  they 
become  even  smaller  in  Parts  1  and  2  for  short  geodesics  and  in 
Part  3  for  large  absolute  latitudes.  It  should  be  noted  that  terms 
containing  powers  of  n  are  in  radians. 

The  accurate  floating  point  calculations  for  short  geodesics 
should  be  applied  not  only  to  the  formulas  of  this  appendix  but,  in 
turn,  also  to  associated  formulas  in  Appendices  A  and  B,  as  illu¬ 
strated  numerically  in  Appendix  D.  The  prescribed  increase  in  deci¬ 
mal  accuracy  in  the  sine  of  a  small  angle,  for  example,  can  be  ob¬ 
tained  not  only  from  the  sine  series,  but  also  from  trigonometric 
tables  by  taking  the  reciprocal  of  the  large  interpolated  cosecant 
of  the  angle.  However,  in  addition  to  sufficient  significant 
digits,  the  table  should  have  intervals  small  enough  for  accurate 
linear  interpolation.  Even  better,  of  course,  is  a  table  of  high 
decimal  accuracy  for  the  small  sines  themselves. 
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APPENDIX  D 


NUMERICAL  ILLUSTRATIONS  OF  INVERSE  AND  DIVtCT 
(Geodesics  of  approx imately  1  end  6,000  miles  for  each) 

The  two  extreme  cest  distances  noted  above  are  chosen  to  illus¬ 
trate  by  calculation  not  or.iy  rhe  basic  computation  forms  of  Appen¬ 
dices  A  ana  B  but  also  the  alternate  formulas  of  Appendix  C.  The 
degree  of  consistency  of  the  answers  has  been  determined  below  by 
checking  each  Inverse  solution  against  the  corresponding  Direct. 

The  resulting  discrepancies,  which  for  each  geodesic  arc  summarized 
at  the  end  of  this  appendix,  therefore  represent  the  combined  errors 
of  the  Inverse  and  Direct. 


Inverse  Solution 
B- 

h 

a.Q  (meters) 
bQ  (meters) 
f 
n 
L 

tan  3, 
tar.  32 

COS 

cos 
p,-r.  gx 
sin  S2 


Long  Geodesic 
+2Cf 
t+ 

+45° 

+106° 

6  378  388.000 
6  356  911.946 
.00  33670  03367 

+106 ' 

.36274  47453 
.99663  29966 
.94006  23275 
.70829  81969 
.34100  26695 
.70591  33545 


Short  Geodesic. 
+45° 

+12° 11’  18" 
+45° 00'36".5 
+12° 12’ 09”. 5 
6  378  388.000 
6  j.'6  911.946 
00  33670  0336/ 
00  16863  40641 
51". 5 

.99663  29966 
.°d698  57825 
.70829  81969 
.708.:  32700 
.70591  33545 
.70603  3681/ 
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Inverse  Solution 

Long  Geodesic 

Short  Geodesic 

a 

.24071  83383 

.49840  21342 

b 

.66584  44515 

.50159  78502 

sin  L 

.96126  16959 

.000  24967  90432 

cos  L 

-.27563  ^3558 

.99999  99688 

COS  0 

.05718  67343 

.99999  99687 

sin  (Bg  -  B1 ) 

.000  17695  69927 

(3g  -  p3  )  radians 

.000  17695  60928 

sin  (0a  “  Pi ) 

.000  17695  60919 

sin2  - 
2 

.000  00001  55849 

sin  0 

.99836  34996 

.000  25016  57049 

0  (radians) 

1.51357  83766 

.000  25016  57075 

c 

.64109  99269 

.50062  20631 

m 

.58899  08837 

.74937  75499 

(S  +  bj  = 

+1.51869  17590 
+  .00080  87665 

-  .00156  22320 

-  .00000  00188 
+  .00000  01279 

+.000  25101  08523 
+.000  00042  05152 
-.000  00063  22699 
-.000  00000  03522 
-.000  00000  07962 

+  .00000  18468 

+.000  00000  10592 

S  (meters) 

9  649  412.505 

1  594.307  213 

P 

i 

t-1 

-i- 

o 

II 

”  +  .00511  33825 

-  .00000  76243 

-  .00001  16616 

”+.000  00084  51448 
-.000  00000  21203 

.000  00000  ooooo 

X  (radians) 

1.85331  48325 

.000  25010  10825 

sin  A. 

.96035  63877 

.000  25010  10799 

cos  X 

-.27877  51927 

.99999  99687 

.  2  X 

Sin  2 

----- 

.000  00001  56376 
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Inverse  Solution 

Long  Geodesic 

Short  Geodesic 

cot  o1!  _2 

1.07455  96453 

.99919  16383 

cot  a3_x 

-.47245  22960 

.99883  88553 

°h  -3 

42°56'30". 03503 

45°01 ’23". 40210 

«a-i 

295° 17 ' 18". 59981 

225°01 ' 59". 82121 

Direct  Check 

Long  Geodesic 

Short  Geodesic 

4-20° 

445° 

h. 

0° 

+12° 11 1 18" 

«i.a 

42°  56 ' 30". 03503 

45°01'23". 40210 

S  (meters) 

9  649  412.505 

1  594.307  213 

aQ  (meters) 

6  378  388.000 

6  378  388.000 

bQ  (meters) 

6  356  911.946 

6  356  911.946 

f 

.00  33670  03367 

.00  33670  03367 

e’2 

.00  67681  70197 

.00  67681  70197 

n 

.00  16863  40641 

tan 

.36274  47453 

.99663  29966 

cos  px 

.94006  23275 

.70829  81969 

sin  g1 

o /.  inn  Q££nc; 

•  J4J.U9  4UU7J 

.70591  33545 

sin  ax_2 

.68125  35334 

.70739  26381 

cos  O'!  _a 

.73204  75552 

.70682  08083 

cos  p0 

.64042  07822 

.50104  49301 

g 

.68817  03286 

.50063  99041 

“a 

.59009  33386 

.25146  93699 

0g  (radians) 

1.51794  02494 

.000  25079  90085 
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Direct  Check 


sin  0£ 


cos  0 


0O  (radians)  = 


sin  0O 
COS  0g 
cot  ff3_x 

ff2-l 
COt  \  = 


X 


(L  -  x)  T  COS  g0  = 


L  (radians) 


L, 

<e 

sin  g2 


cos  p2 
tan  f}2 


tan  B 


Long  Geodesic 
.99860  34425 


Short  Geodesic 


.000  25079  90059 


,05283  14696 


.99999  99685 


.24057  82171 


,49924  27565 


+1.51794  02494 

-  .00081  30002 

-  .00146  29306 
+  .00000  00374 
+  .00000  39825 
+  .00000  25543 


+1.51567  09428 


+.000  25079  90085 
-.000  00042  37199 
.000  00000  00000 
+.000  00000  17897 
.000  00000  000C0 
.000  00000  00000 


+.000  25037  70783 


,99848  09807 


.000  25037  70757 


.05509  74693 


,99999  99687 


-.47245  22450 
295° 17 ’18". 59121 


.99883  88542 
225°01'59". 82132 


-.29028  29979  tan  \  =  .000  25010  10884 


106°11'13". 61256 


51". 58705  146 


-.00511  09099 
+.00000  40853 
+.00000  73513 


-.000  00084  44411 
+.000  00000  21292 
.000  00000  00000 


1.85004  89647 
105°59'59". 99117 
.70591  33687 


.000  24967  90471 
12° 12 '09". 50000  027 


70(390  Q1.Q9Q 

;  /  v/  \j  x-  j  \j  x.  yj  <u  y 


ns 

H 

(D 


99663  30364 


1.00000  00399 
45° OO'OO". 00411 


H* 

9 

<  H* 
ftj  p 
Pi 

c  t-t 
re  ^ 
u 


.70603  86812 
.70817  32700 
.99698  57817 
1.00035  39769 
45°00 ' 36". 49992 
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Direct  Check  Long  Geodesic  Short  Geodesic 

sin3  X  -  -  .000  00001  56376 

2 

sin  (ft,  -  fr)  .  -000  17695  60922 

(Pa  -  px )  radians  -  -  -  -  -  .000  17695  60931 

cos  (Bg  +  Bx )  _____  -.00017  69566 

cos  (Ba  -  Bj)  _____  .99999  99841 

Bg  (improved  value)  _____  45°00 ' 36". 50000  005 


7 


Discrepancies  between 


Inverse  and  Direct 

Long  Geode sit 

:  Short  Geodesic 

ARg 

0". 00411 

0". 00000  005 

AI* 

0". 00883 

0". 00000  027 

■ — ! 

i 

0". 00860 

0". 00011 

In  addition,  the  preceding  Inverse  and  direct  illustrative  examples 
contain  several  common  intermediate  and  secondary  components  whose 
values  can  be  compared.  Also,  since  the  solutions  of  the  long 
geodesic  are  illustrated  by  the  same  numerical  problem  that  was  used 
in  reference  [1]  for  the  earlier  form  of  the  Inverse  method,  oppor¬ 
tunities  for  other  comparisons  are  available.  It  is  apparent  that 
the  extremely  high  positional  accuracy  for  the  short  geodesic  is  due 
to  the  use  of  alternate  formulas  given  in  Appendix  C.  The  azimuth 
error  is  consistent  with  this  positional  error,  in  view  of  the  line's 
shortness.  Comparable  accuracies  are  ajlso  obtainable  at  large 
absolute  latitudes,  but  only  if  interpreted  relative  to  the  in¬ 
creasing  convergence  and  closeness  of  the  meridians  in  polar 
regions . 
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APPENDIX  E 


THEORETICAL  FORMULAS  FOR  HIGHER  ACCURACY 


The  results  of  the  illustrative  numerical  examples  given  in 
Appendix  D  indicate  that  the  formulas  in  Appendices  A  through  C 
provide  sufficient  practical  accuracy.  For  theoretical  purposes, 
however,  the  formulas  could  be  extended  through  f3  and  e'9  terms  or 
beyond.  The  outer  coefficients  of  the  formula  for  (S  4  bQ)  in 
Appendix  A  would  then  include,  for  example,  the  higher  order  combi¬ 
nations  a3,  m3 ,  a3m,  and  am2.  Similar  orderly  extensions  should  be 
expected  for  the  (\  -  L)  4  c  formula  in  Appendix  A  and  the  0Q  and 
(L  -  \)  4-  cos  3q  in  Appendix  B,  except  that  in  the  case  of  the 
latter  two  their  outer  coefficients  will  bear  the  subscript  1,  and 
their  components  and  11^  would  have  to  be  properly  defined  to 
higher  powers  of  e'2.  If  necessary,  appropriate  formulas  in  Appendix 
C  can  also  be  extended. 

In  the  present  appendix,  only  the  (\  -  L)  4-  c  power  series  of 
Appendix  A  will  be  given  to  the  next  higher  order  terms,  since  it 
provides  a  non-iterative  rigorous  solution  for  the  quantity  \  which 
is  required  in  most  of  the  classical  methods  for  calculating  the. 
Inverse  of  long  geodesics.  The  unique  form  of  the  extended  (\  -  L) 

4-  c  power  series  given  below  has  been  derived  from  the  top  of  page 
18  of  [1],  by  substitution  in  terms  of  a,  m,  0,  and  f.  The  series 
is  followed  by  accurate  Inverse  distance  and  azimuth  formulas  taken 
in  large  part  from  pages  14  and  15  of  reference  [1].  The  resulting 
method  of  solution  can  be  used  for  precise  computation  of  Inverse 
problems,  especially  as  a  theoretical  check  on  Direct  or  other  In¬ 
verse  formulas. 

C(f  +  fS  +  f3)  0] 

f3 

f  a  [  -  (j-  4  f3)  sin  0  -  (f3  4  4f3)  02  esc  0 
3f3 

+  (~“)  03  CSC  0  cot  0] 

-  5  f 3  f 8  f 3 

4  m  [_  -  (— —  4-  3f3)  0  4  (“  4  2~ )  sin  0  cos  0 

4  (f3  4  4f3)  0s  cot  0  -  (-^r)  03  esc3  0  -  (f3)  03  cot2  0] 

f3 

4  a3  [(f3)  0  4  (j~)  sin  0  cos  0  4  (f3)  03  esc3  0] 
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-  1 1  c3  n  c3  j:3 

+  rns  r (~15~)  0  '  sin  0  cos  0  +  (j~)  0  cos2  0 

-  (~4~)  02  co,:  0  +  (-— )  sin  0  cos3  0 

2  8 

+  (■—)  03  CSC13  0  +  (2f3)  03  cot2  0] 


where  the  component  quantities  are  again  defined  in  Appendix  A, 
while  some  alternate  definitions  are  found  in  Appendix  C. 

Next,  0Q  is  obtained  in  the  same  manner  as  0,  except  that  the 
value  of  \  obtained  from  above  is  now  to  be  used  in  place,  of  L. 
Then  continue  as  follows: 

cos  =  (b  sin  \)  -  sin  0O  ;  cos  2 a  =  (2a  t  sina  |30)  -  cos  0O« 


-  1  +  T  sin2  Po  -  TT  sin4  +  §T  sin6 

-  ”  sin2  30  -  s in4  j30  +  ^  sin3  0O 


e'4  •  4  a  3e ..e  Q 

1*28  Sin  '  512  S1J  ^ 


D'J  1536  SLn8 


S  =  bQ  (AQ  0O  +  B0  sin  0Q  cos  2cj  -  CQ  sin  20Q  cos  4a 
+  Dn  sin  30o  cos  6a) 

To  complement  the  above  geodetic  distance,  S,  the  azimuths 
O',  „  and  a0  ,  are  obtained  from  formulas  given  in  Appendix  A  or  C. 
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APPENDIX  F 

INTER-RELATIONS  OF  THE  TERMS  OF  THE  POWER  SERIES 


As  noted  earlier,  the  coefficients  a  and  m  in  the  (S  r  b  )  and 
(\  -  L)  -f  c  Inverse  power  series  in  Appendices  A  and  E  display  a 
unique  set  of  product  combinations.  The  identical  simple  pattern 
is  also  repeated  in  the  two  Direct  power  series  in  Appendix  B,  ex¬ 
cept  that  it  occurs  instead  with  the  subscripted  and  .  Al¬ 
though  not  shown  in  this  paper,  even  the  higher  degree  combinations 
(such  as  a2m,  m2a,  a3,  and  m3)  appear  to  enter  in  orderly  fashion 
in  the  further  extension  of  the  power  series.  It  is  of  significant 
importance  that  the  a  and  m  (or  a^^  and  )  combinations  are  com¬ 
pletely  factorable  from  the  power  series  terms,  since  this  permits 
the  latter  to  be  tabulated  as  a  function  of  only  the  variable  0  or 
0g  and  the  parameter  f  or  e's.  Electronic  computer  programming  and 
calculations  also  become  simpler,  whether  for  producing  just  the 
table  or  for  calculating  the  entire  Inverse  or  Direct. 


Another  interesting  inter-relation  of  the  terms  of  the  series 
concerns  the  numerical  coefficients  of  the  powers  of  f  and  e'2. 

It  should  be  noted,  for  example,  that  in  Appendix  B  the  numerical 
coefficients  related  to  the  m^  terms  of  the  0Q  power  series  are: 

111315 
64’  64’  8’  32* 


The  total  of  the  above  four  numbers  is  found  to  be  exactly  zero. 

Upon  closer  inspection,  it  is  found  from  the  power  series  in  Appendi¬ 
ces  A,  B,  and  E  that  the  zero  sum  occurs  with  all  sets  of  terms 
having  m  or  as  one  of  the  factors,  even  for  the  (S  4-  bQ)  series 
in  .Appendix  A,  if  it  is  modified  as  shown  later.  When  different 
powers  of  f  are  present,  the  sum  is  zero  separately  for  the  numeri¬ 
cal  coefficients  of  the  f  terms,  f2  terms,  and  so  forth,  such  as 
in  the  (\  -  L)  4  c  series  in  Appendix  E.  In  all  instances  described, 
the  sum  is  zero  by  virtue  of  the  fact  that  each  term--which  is  a 
function  of  0  (or  0g)--is  first  put  into  a  form  which  satisfies  the 
following  condition:  The  algebraic  sum  of  the  exponents  of  0  and 
sin  0  (after  all  trigonometric  functions  of  0  are  converted  to  sines 
and  cosines)  is  unity.  Actually,  the  above  condition  can  be  (and 
has  been)  satisfied  even  for  the  non-m  and  the  non-m^  series  terms. 
For  very  short  geodesics  (which  of  course  have  a  small  arc  value  0 
and,  therefore,  sin  0  approaches  0  and  cos  0  approaches  unity),  the 
resulting  unity  exponent  implies  that  every  term  is  of  the  small 
order  of  0,  times  its  numerical  coefficient  and  the  proper  power  of 
f  or  e'2.  Since  even  the  omitted  terms  of  the  series  contain  that 
small  order  of  0  (or  0g) ,  the  power  series  converge  to  a  greater 
number  of  decimals  for  short  geodesics.  This  is  shown  by  the  much 
better  positional  consistency  obtained  from  the  numerical  example 
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for  the  short  geodesic  in  Appendix  D.  For  terms  which  have  in  or  nr 
as  one  of  the  coefficients,  the  convergency  for  short  geodesics  is 
even  greater  because  (as  noted  above)  the  sum  of  the  numerical  coef¬ 
ficients  is  zero  separately  fo’*  each  power  of  f  or  e'2,  and  0  or  <p„ 
is  practically  a  common  factor. 

As  for  the  (S  t  b0)  series  mentioned  in  the  preceding  para¬ 
graph,  the  expression  given  in  Appendix  A  can  be  reduced  to  the 
following  form: 

4=  [(l  +  f+f2)0] 


JT<3  _ 

+  (m  cos  0  -  a)  [-  (f  +  f2)  sin  0  +  (_)  02  esc  0] 


f  4.  f2  f  I  f3 

[-  ( — 2 - ^  0  +  ( — 2 - S>  siR  ^  CCS  f/>^ 


fa 

+  (m  cos  0  -  a)2  [-  (— )  sin  0  cos  0J 

+  m2  [ (j^)  0  +  (j^-)  sin  0  cos  0  -  (-— )  sin  0  cos3  0] 

ji2  £2 

+  m(m  cos  0  -  a)  [(-— )  sin  0  cos2  0  -  (y- )  02  esc  0] 

The  compound  coefficient  (m  cos  0  -  a)  is  an  expression  which 
appeared  extensively  in  the  course  of  the  original  derivation  of 
the  Inverse  solution.  As  used  above,  it  causes  the  numerical  coef¬ 
ficients  of  the  terms  with  the  factor  m  to  add  to  zero,  just  like 
the  other  power  series.  It  is  interesting  to  note  that  the  next 
higher  order  extension  of  (S  t  bQ)  continues  to  give  the  proper 
zero  sum  for  the  numerical  coefficients  of  applicable  terms,  when 
the  additional  prescribed  product  combinations  of  the  same 
(m  cos  0  -  a)  and  m  are  used. 


In  conclusion,  it  is  worth  noting  that,  of  the  four  main  power 
series  given  in  Appendices  A  and  B,  only  (S  -f  bQ)  does  not  lend  it¬ 
self  to  completely  factoring  out  the  ellipsoidal  parameter  from  each 
series  of  terms.  The  capability  of  factoring  for  all  four  power 
series  (at  least  to  the  extent  of  the  number  of  terms  given)  may  be 
important.  It  would  mean,  for  example,  that  the  total  value  of 
each  series  of  terms  could  be  tabulated  independently  of  any  specific 
spheroid  flattening  or  eccentricity.  (Of  course,  the  parameters 
would  then  be  made  a  part  of  the  external  coefficients  instead.) 

In  the  (S  -f  bQ)  formula  given  in  the  present  ippendix,  only  the 
terms  whose  coefficient  is  (m  cos  0  -  a)  do  not  lend  themselves  to 
factoring  out  the  function  of  flattening.  Those  terms,  howe’er, 
can  be  represented  as  in  the  following: 
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[  (in  cos  0  -  a)  (1  -  —LZr-r)  ]  [  -  (f  +  f2)  sin  0  ] 

2  sin  0 

where  the  unwanted  portion  of  flattening  has  been  transferred  to 
the  external  coefficient.  This  new  compound  coefficient  may  be  used 
in  place  of  the  previous  (m  cos  0  -  a)  throughout  the  (S  t  b0) 
expression  for  consistency,  since  the  extraneous  f3  terms  which  are 
introduced  are  negligible. 
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APPENDIX  G 

MERIDIONAL  ARC  AS  SPECIAL  CASE  OF  NON-ITERATIVE  INVERSE  AND  DIRECT 


An  interesting  indication  of  the  simplicity  and  rapid  con¬ 
vergence  of  the  non-iterative  inverse  is  to  reduce  it  to  the  special 
case  of  meridional  arc  distances,  for  northern  latitudes  up  to  90° 
from  the  equator.  Since  ^  and  L  are  then  0°,  the  following  result: 


a  =  0,  m  =  1,  0  =  (32  radians. 


Therefore,  such  meridional  distances,  S  ,  become: 

f  0  f2  f  -I  r2 

”  ^0  [(  1  +  2  +  16  )  fts  ■  (  2  +  fd”  ^  sin  3s  cos  3e 


-  (  ~  )  sin  32  cos3  32  ] 


Similarly,  32  can  be  derived  for  the  corresponding  S  by 
letting  p1  and  _2  equal  05  in  the  Direct  solution,  whence: 


ai  =  0,  mj  =  1,  32  -  0O. 

Then  by  substitution  into  the  0Q  power  series,  there  results: 


„  i2  1  1  p  14  p  (S 


i  >4 

— — —  )  sin  0  cos  0 
64  rl 


e'4  £i4 

--7T-  )  sin  0  cos3  0  -  (  —  —  )  0  cos3  0  radians, 
32  HM  8  *M 


where  0  =  (  S  t  b.  )  radians. 

M  M  0 

As  a  check  (the  complete  details  of  which  need  not  be  shown), 
the  above  Inverse  and  Direct  meridional  arc  solutions  were  compared 
mathematically  and  found  to  be  fully  consistent  with  each  other. 
Essentially,  dividing  the  Inverse  meridional  formula  by  b0  produced 
0  as  a  function  of  g2 ,  from  which  sin  0^  and  cos  0^  were  then  ob¬ 
tained  by  expanding  in  series  around  sin  32  and  cos  |32  ,  respectively. 
Substitution  into  the  Direct  meridional  formula  finally  made  the 
right  side  identically  equal  to  the  left  side's  32  ,  up  through  all 
e'4  terms. 
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SUMMARY 


Improved  practical  and  theoretical  formulas  are  presented  for 
the  calculation  of  geodetic  distances,  a” Snaths,  nH  positions  <?n 
spheroid.  The  formulas  are  designed  for  use  with  either  electron! 
Computers  or  desk  calculators.  For  the  latter,  the  formulas  lend 
themselves  to  the  construction  of  useful  interpolation  tables 

The  report  includes  convenient  computation  forms  and  auxiliar 
equal  ions  which  assure  r  high  degree  of  accuracy  for  any  geodetic 
line,  no  matter  ho*  short  or  long,  (up  to  half  or  fully  around  the 
earth)  and  regardless  of  its  orientation  or  location.  Numerical 
examples  Illustrate  the  complete  calculation  procedure. 
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